{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This tutorial will setup the necessary system for running a simple graphene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of manually defining the graphene system with associated atomic coordinates and lattice vectors we use the build-in `sisl` capability of defining the graphene structure with a default atomic distance of $d = 1.42\\,A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene()\n",
    "print(graphene)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some basic information of the geometry is shown above:\n",
    "\n",
    "1. It contains 2 atoms, `na` is *number of atoms*\n",
    "2. It contains 2 orbitals, `no` is *number of orbitals*\n",
    "3. Then there is a list of unique associated atoms. In this case there is 1 unique species, `species: 1` while the subsequent lines specifies the number of times it is occuring in the `Geometry`. \n",
    "4. The last item, `maxR` is the maximal orbital range in the geometry. In this case it is equal to `maxR` as listed for the atomic Carbon species. However, in case there are more than one specie it will be the maximum for the individual species.\n",
    "5. `nsc` specifies the *number of supercells* which is a practical way of dealing with periodic structures.\n",
    "In this case we have periodicity along the first and second lattice vector. And there are 3 images along each of these directions. The primary, and one at $+$ and $-$. It should be apparent that the values in `nsc` are **always** uneven numbers. Note that for molecules it is always `[1, 1, 1]`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to proceed with creating the Hamiltonian for the graphene lattice. But before doing so you should now consider which type of system we are going towards? Are we going to implement a nearest, next nearest or third nearest neighbour interaction?  \n",
    "*HINT*: this should be apparent from the `maxR` variable above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = sisl.Hamiltonian(graphene)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `H` is a Hamiltonian object. It works *equivalently* to a matrix and one may assign elements and extract elements as though it were a matrix, we will return to the intricate details of the Hamiltonian object later.\n",
    "\n",
    "There are only 3 more things that specifies the Hamiltonian. \n",
    "1. `non-zero` is the current number of non-zero elements specified in the Hamiltonian\n",
    "2. `orthogonal` specifies whether an overlap matrix is associated with the Hamiltonian, i.e. `False` for a non-orthogonal basis-set.\n",
    "3. The number of spin-components are specified via the `Spin` class which in this case is an unpolarized system. I.e. only one spin-component.\n",
    "4. The `Geometry` class is repeated because the Hamiltonian is associated with the `Geometry` passed when instantiating the Hamiltonian (`sisl.Hamiltonian(graphene)`).\n",
    "\n",
    "You are now ready to add matrix elements to the Hamiltonian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explicitly specifying the matrix elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this small sub-section we will specify the required matrix elements manually. You should already now know that we only plan on adding nearest neighbour interactions.\n",
    "\n",
    "We will proceed with an on-site of $0\\,\\mathrm{eV}$ and a coupling element of $t=-2.7\\,\\mathrm{eV}$.  \n",
    "First we set the on-site elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H[0, 0] = 0.0\n",
    "H[1, 1] = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to set the coupling elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H[0, 1] = -2.7\n",
    "H[1, 0] = -2.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will only couple the first and second atom in the primary unit-cell. But we also require couplings from the primary unit-cell to the neighbouring supercells. Remember that `nsc = [3, 3, 1]`.  \n",
    "In this case we know that the first atom couples to the `(-1, 0)` and `(0, -1)`, while the second atom couples to `(1, 0)` and `(0, 1)`.  \n",
    "The coupling elements may then be specified via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H[0, 1, (-1, 0)] = -2.7\n",
    "H[0, 1, (0, -1)] = -2.7\n",
    "H[1, 0, (1, 0)] = -2.7\n",
    "H[1, 0, (0, 1)] = -2.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all matrix elements are set, i.e. 2 on-site and 6 nearest neighbour couplings, lets assert this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We find 8 non-zero elements, as there should be. Remark that even though we set the on-site terms to $0$, they are interpreted as non-zero elements due to explicitly setting them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After having setup the Hamiltonian, we may easily calculate the eigenvalues at any $\\mathbf k$ (in reduced coordinates $\\mathbf k\\in]-0.5:0.5]$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Gamma:\", H.eigh())\n",
    "print(\"K:\", H.eigh(k=[2./3,1./3,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looping the atoms and orbitals in the Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above specification of the elements is tedious for anything with more elements than 2 orbitals. There are easier and better ways to set all the coupling elements.  \n",
    "In the following we will introduce a new function:\n",
    " \n",
    "     Geometry.close\n",
    "\n",
    "which is a *very* convenient function to return all atomic indices within a certain radii from a given coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ia, io in H.geometry.iter_orbitals(local=False):\n",
    "    # This loops over all atoms and the orbitals\n",
    "    # corresponding to the atom.\n",
    "    # In this case the geometry has one orbital per atom, hence\n",
    "    #   ia == io\n",
    "    # in all cases.\n",
    "    \n",
    "    # In order to figure out which atoms atom `ia` is connected\n",
    "    # to, we must find those atoms.\n",
    "    # To do this we access the geometry attached to the \n",
    "    # Hamiltonian (H.geom)\n",
    "    # and use a function called `close` which returns ALL \n",
    "    # atomic indices within certain ranges of a given point or atom\n",
    "    idx = H.geometry.close(ia, R = [0.1, 1.43])\n",
    "    # the argument R has two entries:\n",
    "    #   0.1 and 1.43\n",
    "    # Each value represents a radii of a sphere.\n",
    "    # The `close` function will then return\n",
    "    # a list of equal length of the R argument (i.e. a list with\n",
    "    # two elements).\n",
    "    # idx[0] is the first element and is also a list\n",
    "    # of all atoms within a sphere of 0.1 AA of atom `ia`.\n",
    "    # This should obviously only contain the atom it-self.\n",
    "    # The second element, idx[1], contains all atoms within a sphere\n",
    "    # with radius of 1.43 AA, but not including those within 0.1 AA.\n",
    "    # In this case this is then all atoms that are the nearest neighbour\n",
    "    # atoms.\n",
    "   \n",
    "    # Now we know the on-site atoms (idx[0]) and the nearest neighbour\n",
    "    # atoms (idx[1]), all we need to do is to set the Hamiltonian\n",
    "    # elements:\n",
    "\n",
    "    # on-site (0. eV)\n",
    "    H[io, idx[0]] = 0.\n",
    "   \n",
    "    # nearest-neighbour (-2.7 eV)\n",
    "    H[io, idx[1]] = -2.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above loop is equivalent to the previously explicitly set values, so printing the structure will yield the same information, we have just specified all values again.  \n",
    "The above loop is also *only* 4 lines of code, irrespective of the number of atoms in the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After having setup the Hamiltonian, we may easily calculate the eigenvalues at any $\\mathbf k$ (in reduced coordinates $\\mathbf k\\in]-0.5:0.5]$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Gamma:\", H.eigh())\n",
    "print(\"K:\", H.eigh(k=[2./3,1./3,0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may also create a band-structure of the Hamiltonian.  \n",
    "First we define the path we wish to take in the Brillouin zone. Start at $\\Gamma$, continue to $K$, then $M$ and finish where we started.  \n",
    "We will calculate the band-structure with a total of `301` points and we explicitly name the special points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band = sisl.BandStructure(H, [[0., 0.], [2./3, 1./3],\n",
    "                              [1./2, 1./2], [1., 1.]], 301,\n",
    "                              [r'$\\Gamma$', 'K', 'M', r'$\\Gamma$'])\n",
    "eigs = band.apply.array.eigh()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now `eigs` contains all the eigenvalues of the Hamiltonian object for all the $k$-points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the tick-marks and the linear k points\n",
    "xtick, xtick_label = band.lineartick()\n",
    "lk = band.lineark()\n",
    "plt.plot(lk, eigs)\n",
    "\n",
    "plt.ylabel('Eigenspectrum [eV]')\n",
    "plt.gca().xaxis.set_ticks(xtick)\n",
    "plt.gca().set_xticklabels(xtick_label)\n",
    "\n",
    "# Also plot x-major lines at the ticks\n",
    "ymin, ymax = plt.gca().get_ylim()\n",
    "for tick in xtick:\n",
    "    plt.plot([tick,tick], [ymin,ymax], 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "- Go back to `In [4]` and define a potential difference between the two atoms (sub-lattices, difference should be more than `0.5 eV`).\n",
    "   What happens to the band structure?\n",
    "- Change the hopping element.\n",
    "   What happens to the band-structure?\n",
    "- Learn how to use `Geometry.close`. Determine the number of indices returned by the 4 below calls, and why? It may help to draw the graphene lattice and geometry on a piece of paper.\n",
    "\n",
    "    1. Do 3 `close` calls at coordinate $(1, 0, 0)$ for 3 different radius, $R \\in \\{0.5, 1.45, 1.7\\}$ (`R = 0.5`...).  \n",
    "    2. Do 1 `close` call at coordinate $(1, 0, 0)$ for 3 different radius, `R = [0.5, 1.45, 1.7]`.\n",
    "         Note, that in this case a list of 3 lists are returned (one per radius)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## Learned lessons\n",
    "\n",
    "- Usage of `sisl.geom` class to create default geometries\n",
    "- Creation of a Hamiltonian object (`Hamiltonian`)\n",
    "- Assigning elements to the Hamiltonian matrix\n",
    "- Retrieving atoms within spherical radius of points or atomic indices (`Geometry.close`)\n",
    "- Creating a band-structure (`BandStructure`)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
